#include "feElement.h"
#include "feSmartPtr.h"
#include "feMvOps.h"
#include "feMatrix.h"
#include "feVector.h"
#include "feMaterial.h"

double feLinearElement1d::WeightRightCoeff = 1;
double feLinearElement1d::StiffRightCoeff  = 0;




feQuadraticElement1d::feQuadraticElement1d(int iStartNode, feMesh* iMesh, feSharedPtr<feMaterial> iMat)
{
	StartNode = iStartNode;
	Mesh = iMesh;
	Material = iMat;
	Lookup1 = -1; 
	Lookup2 = -1;
	Lookup3 = -1; 

	WeightRightCoeff = 1;
	StiffRightCoeff  = 0;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
fePoint feQuadraticElement1d::NodeCoords(int iLocalN)
{
	if(iLocalN != 1)
		return Mesh->GetNode(NodeIndex(iLocalN));

	return (Mesh->GetNode(NodeIndex(0)) + Mesh->GetNode(NodeIndex(2))) / 2;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
int feQuadraticElement1d::NodeIndex(int iLocalN) 
{
	return StartNode+iLocalN;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feQuadraticElement1d::GetStiffness(feDenseMatrix **oMatrix)
{
	//calc h
	fePoint p1 = Mesh->GetNode(StartNode);
	fePoint p2 = Mesh->GetNode(StartNode + 2);
	double h = p2.X - p1.X;


	//create matrix
	double m[3][3] = 
	{
		{7/(3*h), -8/(3*h), 1/(3*h) },
		{-8/(3*h), 16/(3*h), -8/(3*h) },
		{1/(3*h), -8/(3*h), 7/(3*h) }
	};

	
	//allocate if necessary
	if(!(*oMatrix))
		(*oMatrix) = new feDenseMatrix(3);
	else if((*oMatrix)->GetSize() != 3)
	{
		delete (*oMatrix);
		(*oMatrix) = new feDenseMatrix(3);
	}


	//fill matrix
	double *pm[3];
	pm[0] = (double *)&m[0];
	pm[1] = (double *)&m[1];	
	pm[2] = (double *)&m[2];
	(*oMatrix)->Fill((double**)pm);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feQuadraticElement1d::GetWeight(feDenseMatrix **oMatrix)
{
	//calc h
	fePoint p1 = Mesh->GetNode(StartNode);
	fePoint p2 = Mesh->GetNode(StartNode + 2);
	double h = p2.X - p1.X;


	//create matrix
	double m[3][3] = 
	{
		{4*h/30, 2*h/30, -1*h/30 },
		{2*h/30, 16*h/30, 2*h/30 },
		{-1*h/30, 2*h/30, 4*h/30 }
	};


	//allocate if necessary
	if(!(*oMatrix))
		(*oMatrix) = new feDenseMatrix(3);
	else if((*oMatrix)->GetSize() != 3)
	{
		delete (*oMatrix);
		(*oMatrix) = new feDenseMatrix(3);
	}


	//fill matrix
	double *pm[3];
	pm[0] = (double *)&m[0];
	pm[1] = (double *)&m[1];	
	pm[2] = (double *)&m[2];	
	(*oMatrix)->Fill((double**)pm);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feQuadraticElement1d::GetRight(feVector *oVector, int iLayer)
{
	//calc h
	fePoint p1 = Mesh->GetNode(StartNode);
	fePoint p11 = Mesh->GetNode(StartNode+1);
	fePoint p2 = Mesh->GetNode(StartNode + 2);
	double h = p2.X - p1.X;


	//calc f[i]
	double f1 = Material->GetF(&p1, iLayer, &Lookup1);
	double f2 = Material->GetF(&p11, iLayer, &Lookup2);
	double f3 = Material->GetF(&p2, iLayer, &Lookup3);


	//allocate if necessary
	if(oVector->GetSize() != 3)
	{
		oVector->Free();
		oVector->Allocate(3);
	}


	double s = StiffRightCoeff;
	double w = WeightRightCoeff;
	//fill vector
	oVector->Data[0] = s*(7*f1 - 8*f2 + 1*f3)/(3*h) +  w*(4*f1 + 2*f2 - f3)*h/30;
	oVector->Data[1] = s*(-8*f1 + 16*f2 - 8*f3)/(3*h) +  w*(2*f1 + 16*f2 + 2*f3)*h/30;
	oVector->Data[2] = s*(1*f1 - 8*f2 + 7*f3)/(3*h) +  w*(-f1  + 2*f2 + 4*f3)*h/30;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feQuadraticElement1d::GetIndeces(int *oIndeces)
{
	//fill vector
	oIndeces[0] = StartNode;
	oIndeces[1] = StartNode+1;
	oIndeces[2] = StartNode+2;
}



//---------------------------------------------------------------------------------------------------------------------------------------------
//feLinearElement1d
//---------------------------------------------------------------------------------------------------------------------------------------------
feLinearElement1d::feLinearElement1d(int iStartNode, feMesh* iMesh, feSharedPtr<feMaterial> iMat)
{
	StartNode = iStartNode;
	Mesh = iMesh;
	Material = iMat;
	Lookup1 = -1;
	Lookup2 = -1;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
fePoint feLinearElement1d::NodeCoords(int iLocalN)
{
	return Mesh->GetNode(NodeIndex(iLocalN));
}


//---------------------------------------------------------------------------------------------------------------------------------------------
int feLinearElement1d::NodeIndex(int iLocalN) 
{
	return StartNode+iLocalN;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLinearElement1d::GetStiffness(feDenseMatrix **oMatrix)
{
	//calc h
	fePoint p1 = Mesh->GetNode(StartNode);
	fePoint p2 = Mesh->GetNode(StartNode + 1);
	double h = p2.X - p1.X;


	//create matrix
	double m[2][2] = 
	{
		{1/(h), -1/(h)},
		{-1/(h), 1/(h) }
	};


	//allocate if necessary
	if(!(*oMatrix))
		(*oMatrix) = new feDenseMatrix(2);
	else if((*oMatrix)->GetSize() != 2)
	{
		delete (*oMatrix);
		(*oMatrix) = new feDenseMatrix(2);
	}


	//fill matrix
	double *pm[2];
	pm[0] = (double *)&m[0];
	pm[1] = (double *)&m[1];	
	(*oMatrix)->Fill((double** )pm);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLinearElement1d::GetWeight(feDenseMatrix **oMatrix)
{
	//calc h
	fePoint p1 = Mesh->GetNode(StartNode);
	fePoint p2 = Mesh->GetNode(StartNode + 1);
	double h = p2.X - p1.X;


	//create matrix
	double m[2][2] = 
	{
		{2*h/6, 1*h/6 },
		{1*h/6, 2*h/6 }
	};


	//allocate if necessary
	if(!(*oMatrix))
		(*oMatrix) = new feDenseMatrix(2);
	else if((*oMatrix)->GetSize() != 2)
	{
		delete (*oMatrix);
		(*oMatrix) = new feDenseMatrix(2);
	}


	//fill matrix
	double *pm[2];
	pm[0] = (double *)&m[0];
	pm[1] = (double *)&m[1];	
	(*oMatrix)->Fill((double** )pm);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLinearElement1d::GetRight(feVector *oVector, int iLayer)
{
	//calc h
	fePoint p1 = Mesh->GetNode(StartNode);
	fePoint p2 = Mesh->GetNode(StartNode + 1);
	double h = p2.X - p1.X;


	//calc f[i]
	double f1 = Material->GetF(&p1, iLayer, &Lookup1);
	double f2 = Material->GetF(&p2, iLayer, &Lookup2);


	//allocate if necessary
	if(oVector->GetSize() != 2)
	{
		oVector->Free();
		oVector->Allocate(2);
	}


	double s = feLinearElement1d::StiffRightCoeff;
	double w = feLinearElement1d::WeightRightCoeff;
	
	//fill vector
	oVector->Data[0] = s*(f1-f2)/h + w*(2*f1 + f2)*h/6;
	oVector->Data[1] = s*(-f1+f2)/h + w*(1*f1 + 2*f2)*h/6;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
void feLinearElement1d::GetIndeces(int *oIndeces)
{
	//fill vector
	oIndeces[0] = StartNode;
	oIndeces[1] = StartNode+1;
}
